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Abstract 


A  least-squares  mimmization  technique  for  estimating  principal  components  (empirical 
orthogonsd  function  expansion  coefficients)  from  incomplete  ocean  profile  data  is  developed, 
and  its  use  for  interpolating  or  extrapolating  these  profiles  is  described. 
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Estimation  of  EOF  Expansion  Coefficients  from  Incomplete  Data 


1.  Introduction 


Among  the  earliest  published  reports  on  the  use  of  empirical  orthogonal  functions  (EOFs) 
in  oceanography  are  included  such  articles  as  Moore  (1974),  Kundu,  et  al.  (1975)  and  Davis 
(1976).  Since  Aen,  the  use  of  EOFs  has  become  increasingly  common  for  such  purposes 
as  describing  and  quantifying  oceanic  variability  (e.g.,  Halliwell  and  Mooers,  1979),  in 
examining  ocean  dynamics  (e.g.,  Kundu,  et  al.,  1975),  in  providing  compact  characterization 
of  oceanographic  data  (e.g.,  Tolstoy,  et  al.,  1991),  and  in  examining  interrelationships 
between  surface  and  subsurface  parameters  (e.g.,  Carnes,  et  al.,  1990).  In  this  note  we 
present  a  technique  for  using  EOFs  for  interpolating  or  extrapolating  missing  data  in  ocean 
profiles.  This  technique  should  prove  particularly  useful  for  extrapolating  shallow  ocean 
data  profiles  --  say,  from  shallow  CTD  (conductivity-temperature-depth)  casts  or  from 
expendable  probes  -  to  deeper  depths. 

The  theoiy  behind  EOF  computations  is  straightforward.  (See,  for  example,  Lagerloef  and 
Bernstein,  1988,  for  a  particularly  compact  and  lucid  description.)  As  an  aid  to  under¬ 
standing  our  technique  described  in  section  2,  we  present  a  summary  of  the  EOF 
computation  procedures  here. 

Suppose  we  have  N  profiles  with  measurements  at  K  discrete  depths.  Each  profile  may  be 
considered  a  vector  f^,  with  components  f|y,.  A  set  of  perturbation  profiles  f  ^  is  formed  by 
subtracting  the  mean  profile  T  from  each  of  the  observed  profiles.  From  these  profiles  the 
K  X  N  matrix  [F]  is  aeated,  having  columns  which  are  the  N  perturbation  profiles  r„.  Next 
the  K  X  K  covariance  matrix  [R]  is  calculated  from  the  [F]  matrix  accor^g  to 

m  «  (l/N)  [FJ  [Ff, 
where  [F]^  is  the  transpose  of  [F]. 

Then  the  eigenvalues  and  eigenvectors  are  found  for  the  matrix  equation 

[R]  [E]  =  [L]  [E]. 


[L]  is  the  matrix  of  eigenvalues 


IL]  = 


X 
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and  [E]  is  the  matrix  of  eigenvectors 


[E]  =  [e^  62 .  • .  e„  •  • .  CiJ 

where  is  a  column  vector  representing  the  mth  eigenmode.  A  maximum  of  K 
eigenmodes  may  be  calculated. 

The  original  data  may  be  expanded  in  terms  of  these  new  basis  functions  e„  according  to 
the  relationship 

[F]  =  [E]  [C]  (1) 

where  the  c,^  elements  of  [C]  are  known  as  the  expansion  coefficients,  the  modal 
amplitudes,  or  the  principal  components.  The  coefficient  matrix  [C]  is  computed  from 

[C]  =  [Ef  [F],  (2) 

The  trace  of  the  matrix  [R]  is  the  total  variance  of  the  data  set.  Since  Tr  [R]  =  z  a„,  the 
relative  value  of  each  is  the  fraction  of  the  variance  of  the  data  represented  by  the 
associated  eigenvector.  If  the  eigenvalues  are  ordered  by  size  (that  is,  by  fraction  of  the 
variance  explained  by  the  corresponding  eigenvector),  it  is  usually  found  that  only  a  limited 
number  of  eigenmodes  are  necessary  to  account  for  a  very  large  fraction  of  the  variance,  say 
3  to  5  (e.g.,  Kundu,  et  al.,  1975;  Servain  and  Legler,  1986;  Fukumori  and  Wunscb,  1991). 
For  many  oceanographic  applications  it  is  satisfactory  to  retain  only  those  first  M 
eigenmodes  that  account  for,  say,  90  or  95%  of  the  variance  and  subsequent  profiles  (the 
"synthetic  profiles,"  as  they  were  called  in  Carnes,  et  al.,  1990)  are  constructed  using  only 
M  eigenvectors. 

If  a  given  profile  f  has  data  at  each  of  the  K  depths,  estimation  of  the  principal  components 
is  straightforward  and  is  given  by  eq.  (2)  above.  However,  if  the  profile  is  incomplete,  a 
modified  procedure  must  be  used.  In  th^  note  we  present  one  such  procedure  based  upon 
least  squares  minimizatioiL 

2.  Tlieoiy 

According  to  eq.  (1)  above,  a  "synthetic  profile"  p  may  be  constructed  from  M  pre¬ 
determined  eigenvectors  according  to  the  equations 

Pk  ~  ®kiii  <^111 

where  k  is  the  depth  index,  k  =  1,  . . . ,  K.  If  all  depths  are  present,  eq.  (2)  shows  the 
coefficients  Cg,  are  calculated  from  the  relationship 
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Cm  =  p’k* 

In  practice,  however,  some  of  the  K  depth  values  may  be  missing;  the  profile  may  not  start 
close  enough  to  the  surface,  or  it  may  have  terminated  at  too  shallow  a  depth  or  there  may 
have  been  noise  contamination  within  a  certain  depth  range.  If  only  an  occasional  value  is 
missing,  it  may  be  reconstructed  by  interpolation  and  the  expansion  coefficients  may  then 
be  computed  from  a  complete  data  set.  However,  if  interpolation  is  not  satisfactory  or 
feasible,  a  least-squares  minimization  approach  may  be  used  in  which  we  determine  those 
principle  components  c„  which  minimize  the  error  c,  the  sum  of  the  squares  of  the 
differences  between  a  perturbation  profile 

P’k=  Pk  -  Tk 

and  its  "synthetic"  reconstruction: 

«  =  2k  (p’k ' Cm  Ckm)^* 

This  involves  solving  the  system  of  equations 


Ckl  Ckl  Ckl  Ck2  •  •  •  Ckl  ^kM 
Ck2  Ckl 

• 

_ 1 

= 

p’k  Ckl  “ 

P’k  Ck2 
• 

CkM  Ckl . C|^  e^M  ^ 

- 1 

•  ^ 

_ 1 

» 

p’k  CkM_ 

(3) 


The  resulting  set  of  coefficients  {Cg,}  constitute  the  desired  estimates  of  the  true  coefficients. 
3.  Example 

As  an  example  of  the  application  of  this  technique,  consider  the  28  temperature  profiles 
shown  in  Fig.  la.  Table  1  lists  the  measurement  depths  of  the  profiles.  Fig.  lb  shows  the 
mean  temperature  profile.  The  data  represent  a  set  of  CTD  temperature  profiles  taken  in 
July  during  an  experiment  in  the  Sargasso  Sea,  but  any  other  comparable  data  set  would  do 
as  well  for  describing  our  technique.  Shown  in  Fig.  2  are  the  first  five  EOF  modes  (i.e., 
through  65),  which  jointly  explain  98%  of  the  variance  apparent  in  Fig.  Ic.  The  fint  five 
principal  components  for  each  perturbation  profile  in  Fig.  Ic  as  estimated  from  eq.  (2)  are 
given  in  Fig.  3a,  and  the  root-mean-square  (RMS)  difference  between  the  true  profiles  and 
the  synthetic  profiles  generated  with  these  five  principal  components  and  the  five 
eigenfunctions  in  Fig.  2  is  given  in  Fig.  3b.  Using  the  five  eigenfunctions,  the  maximum 
RMS  deviation  is  about  0.25  *C  near  50  m,  which  is  located  somewhat  below  the  base  of 
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the  mixed  layer  at  around  30  m.  Otherwise  the  mean  deviation  is  typically  less  than  0.15 
•C 

First  consider  the  impact  on  the  coefficient  values  of  having  missing  data  near  the  surface. 
We  estimated  the  principal  components  for  six  separate  cases  with  increasingly  large 
amounts  of  missing  data,  as  described  in  Table  2.  The  resulting  coefficients,  normalized  by 
their  true  values  computed  from  the  complete  profiles  (Fig.  3a),  are  shown  in  Fig.  4,  and 
the  RMS  differences  between  the  true  profiles  and  their  corresponding  synthetic  profiles 
constructed  from  five  modes  are  given  in  Fig.  5.  Not  surprisingly,  as  the  range  of  the 
missing  data  points  becomes  larger,  in  general  the  deviation  of  the  estimated  coefficients 
from  their  true  values  becomes  larger,  as  does  the  RMS  difference  across  the  depth  range 
of  the  missing  data.  However,  the  coefficients  of  certain  modes  are  much  less  sensitive  than 
the  coefficients  of  other  modes. 

The  coefficients  of  the  first  mode  are  least  sensitive  to  data  deletion,  but  by  the  time  the 
first  100  m  of  data  are  removed  (Case  5),  the  individual  coefficients  begin  to  deviate 
noticeably  from  their  true  values.  Mode  4  coefficients  were  the  most  sensitive,  followed  by 
mode  S,  then  mode  3,  and  mode  2.  Tliis  correlates  well  in  a  qualitative  sense  with  the 
location  of  extrema  in  the  eigenfunctions  of  Fig.  2.  However,  even  with  the  upper  100  (Case 
5)  to  150  (Case  6)  m  of  data  removed,  the  RMS  deviation  is  less  than  1.7  ®C.  Notice  that 
the  RMS  difference  below  the  depth  of  the  missing  data  actually  is  less  than  it  is  when  using 
the  full  profiles.  This  is  because  when  complete  profiles  are  used,  a  trade-off  must  occur 
in  the  cdculation  of  the  principal  components,  where  the  various  eigenfunctions  must  be  fit 
as  well  as  possible  over  the  fiiU  depth  range.  When  the  highly  structured  upper  part  of  the 
profile  is  missing,  the  remaining  portion  may  be  better  fit  with  the  eigenfunctions. 

Next  consider  what  occurs  if  increasing  amounts  of  data  are  missing  at  the  bottom  of  the 
profiles.  The  six  cases  considered  are  described  in  Table  3,  the  normalized  coefficients  are 
in  Fig.  6,  and  the  RMS  differences  are  in  Fig.  7.  Because  there  is  so  much  less  structure 
at  the  bottom  of  the  profiles,  the  impact  of  missing  data  is  less  deleterious.  In  general,  as 
the  amount  of  missing  data  becomes  larger,  the  deviation  of  the  estimated  coefficients  from 
their  true  values  also  becomes  larger.  However,  the  deviations  are  much  smaller  in 
comparison  with  the  estimated  coefficients  constructed  with  near-surface  data  points 
removed.  RMS  differences  are  also  much  smaller;  less  than  0.3  in  Case  12  when  1400 
m  of  data  were  removed  from  the  bottom  of  each  profile. 

Again  the  coefficients  of  the  first  mode  are  by  far  the  least  sensitive  to  the  deletion  of  data; 
they  do  not  vary  significantly  from  their  true  values  even  when  a  large  number  of  near¬ 
bottom  data  points  are  missing  from  the  profiles.  Modes  2  and  5  are  the  most  sensitive,  but 
their  deviations  from  their  respective  true  values  are  much  less  than  when  the  surface  v^ues 
are  removed.  Deviations  in  mode  3  increase  quite  slowly  as  the  amount  of  missing  data 
increases. 

The  RMS  difference  increases  as  the  profiles  are  progressively  truncated.  The  most 
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sensitive  part  of  the  profiles  lies  between  about  600  to  1200  m,  where  mode  2  in  particular 
has  considerable  structure.  Deletion  of  data  below  this  depth  range  has  only  a  modest 
impact  upon  the  resulting  synthetic  profiles.  The  RMS  difference  above  the  truncation 
depths  again  tends  to  be  somewhat  less  than  the  RMS  difference  computed  using  full  depth 
profiles. 

Finally,  we  deleted  increasingly  larger  amounts  of  data  in  the  depth  ranges  centered  around 
300  m  (a  minimum  in  the  variability  seen  in  Fig.  Ic)  and  around  600  m  (a  maximum  in  the 
variability),  as  described  in  Tables  4  (Cases  13,  14,  and  15)  and  5  (Cases  16,  17,  and  18). 
The  normalized  coefficients  and  RMS  differences  are  given  in  Figs.  8  and  9  and  Figs.  10  and 
11.  Both  of  these  sets  of  results  show  deletion  of  data  increases  the  RMS  differences 
between  the  synthetic  profiles  and  the  data  over  that  computed  with  complete  data; 
however,  the  increase  in  the  RMS  is  only  in  the  vicinity  of  the  missing  data  and  it  is 
remarkably  small.  Case  18  was  computed  with  data  between  500  -  700  m  deleted,  but  the 
RMS  difference  at  500  m  increased  less  than  0.1  *C  to  about  0.2  “C. 

4.  Coefficient  Estimation  Software 

To  implement  the  theory  described  above,  a  FORTRAN  subroutine  ESTCOE  was  written 
and  is  presented  in  Appendix  A.  ESTCOE  requires  the  input  of  two  data  files;  one 
containing  the  profiles  themselves  and  the  other  containing  the  mean  profile  and 
eigenvectors  associated  with  the  area  of  interest.  For  each  profile,  the  system  of  equations 
given  in  eq.  (3)  is  solved  for  the  set  of  unknown  estimated  coefficients  {c„}.  The  synthetic 
profile  is  then  constructed  using  the  estimated  coefficients,  the  input  eigenvectors,  and  the 
mean  profile. 

The  output  coefficients  are  stored  in  a  two-^iimensional  variable,  Coef(MaxProfiles,  Maxz), 
where  the  maximum  number  of  profiles  and  maximum  number  of  depths  have  been  set  (by 
a  parameter  statement)  as  MaxProfiles  =  100  and  Maxz  =  50,  respectively.  The  output 
synthetic  profiles  are  also  stored  in  a  two-dimensional  variable,  pp(MaxProfiles,  Maxz). 

Other  parameters  are  returned  from  ESTCOE.  A  full  description  of  these  parameters  as 
well  as  a  more  complete  description  of  the  algorithm,  format  of  input  files,  output 
parameters,  an  example  driver  program,  and  the  FORTRAN  code  itself  is  found  in 
Appendix  A. 

5.  Concluding  Remarks 

We  have  described  a  least-squares  minimization  technique  for  estimating  principal 
components  fi’om  incomplete  ocean  profile  data.  An  important  application  of  this  technique 
is  to  compute  from  these  estimated  principal  components  a  synthetic  profile  that  interpolates 
or  extrapolates  the  missing  data  values.  This  synthetic  proffie  can  then  be  used  in  place  of 
the  original  incomplete  profile,  or  it  can  be  merged  with  the  original  profile  in  some  fashion 
to  fill  in  the  missing  data.  In  particular,  the  technique  lends  itself  to  the  problem  of 


extrapolating  shallow  ocean  data  profiles  ~  say,  from  shallow  CTD  casts  or  from  expendable 
probes  such  as  XBTs  (expendable  bathythermographs)  or  AXBTs  (the  air-launched  XBT)  - 
-  to  a  deeper  preselected  depth.  A  commonly  used  alternative  is  to  use  climatological 
values  below  the  maximum  data  depth,  with  some  merging  procedure  to  combine  the  two. 
The  technique  presented  here,  however,  makes  use  of  the  information  contained  in  the 
structure  of  the  observed  profile  and  the  correlations  between  that  upper  level  structure  and 
the  deeper  structure  to  produce  what  will,  in  many  cases,  be  a  better  estimate  of  the  deeper 
structure  than  can  be  made  with  climatolo^  alone. 

To  illustrate  this,  we  computed  the  RMS  difference  between  our  28  temperature  profiles  in 
Fig.  la  and  the  appropriate  GDEM  (Teague,  et  al.,  1990)  and  Levitus  (Levitus,  1982) 
climatological  profiles  (Fig.  12).  We  also  simulated  the  use  of  T-7  XBTs  by  cutting  our  28 
profiles  off  below  700  m,  estimating  the  principal  components  as  described  above,  and 
computing  the  synthetic  profiles  from  these  principal  components  (this  is  Case  10  of  Fig.  7). 
The  RMS  difference  for  Case  10  between  the  synthetic  profiles  and  the  data  is  also 
reproduced  on  Fig.  12.  The  RMS  difference  for  the  synthetic  profiles  is  substantially  lower 
than  the  RMS  differences  for  either  of  the  climatologies,  particularly  below  700  m.  Hence, 
merging  the  synthetic  profiles  with  the  700  m  deep  measured  profiles  would  give  a  better 
depth-extended  data  set  than  using  either  climatology.  These  profiles  happened  to  come 
from  the  Sargasso  Sea  region,  but  we  feel  confident  our  results  would  hold  elsewhere  as 
well. 
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Table  1.  Standard  depths  (m)  used  in  the  analysis. 
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Table  2.  Depths  at  which  data  values  were  eliminated  for  Cases  1-6. 

Case  1:  0  m 

Case  2:  0,  10,  20  m 

Case  3:  0,  10,  20,  30  m 

Case  4;  0,  10,  20,  30,  50  m 

Case  5:  0,  10,  20,  30,  50,  75,  100  m 

Case  6:  0,  10,  20,  30,  50,  75,  100,  125,  150  m 


Table  3.  Depths  at  which  data  values  were  eliminated  for  Cases  7  -  12. 

Case  7:  2000,  1750  m 

Case  8:  2000,  1750,  1500,  1400, 1300,  1200,  1100  m 

Case  9:  2000,  1750,  1500,  1400, 1300,  1200,  1100,  1000,  900  m 

Case  10:  2000,  1750,  1500,  1400,  1300,  1200,  1100,  1000,  900,  800  m 

Case  11:  2000,  1750,  1500,  1400,  1300,  1200,  1100,  1000,  900,  800,  700  m 

Case  12:  2000,  1750,  1500,  1400,  1300,  1200,  1100,  1000,  900,  800,  700,  600  m 


Table  4.  Depths  at  which  data  values  were  eliminated  for  Cases  13  -  15. 

Case  13:  300  m 

Case  14:  250,300  m 

Case  15:  200,  250,  300,  400  m 


Table  5.  Depths  at  which  data  values  were  eliminated  for  Cases  16  -  18. 


Case  16:  600  m 
Case  17:  500,  600  m 
Case  18:  500,  600,  700  m 
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Depth  (m) 


28  Temperature  Profiles  Mean  Profile  Perturbation  Profiles 


Temperature  (°C)  Temperature  (°C)  Temperature  (°C) 


Figure  1.  The  28  temperature  profiles  used  in  the  analysis,  a;  overplot  of  all  profiles;  b:  mean 
profile;  c:  oveiplot  of  perturbation  profiles. 
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Depth  (m)  Depth  (m) 


Mode  1  (78.6%)  Mode  2  (11.6%)  Mode  3  (4.5%) 


Mode  4  (1.9%) 


Mode  5  (1.3%) 


Amplitude  Amplitude 


Figure  2.  The  first  five  eigenvectors  computed  fi-om  the  data  in  Fig.  1,  The  number  in 
parentheses  is  the  percent  of  the  total  variance  that  is  described  by  that  particular 
eigenvector. 
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Amplitude 


First  5  principal  components 
for  28  temperature  profiles 


Figure  3.  a:  The  expansion  coefficient  values  versus  mode  number  for  coefficients 
calculated  using  full  depth  profiles;  b:  the  RMS  temperature  difference  over  all  28  profiles 
between  the  real  and  the  synthetic  profiles  computed  with  five  modes. 
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Case  1 


Case  2 


Case  3 


Case  5 


Mode 


Mode 


Figure  4.  Hie  normalized  estimated  coefficient  amplitudes  versus  mode  number  for  Cases 
1  -  6  in  which  data  in  the  upper  portion  of  the  profiles  were  deleted.  Table  2  describes  the 
data  deleted  for  each  case. 
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Depth  (m)  Depth  (m) 


Case  1  Case  2  Case  3 


Case  4 


Case  5 


Case  6 


RMS  Difference  (°C)  RMS  Difference  (°C)  RMS  Difference  (^C) 


Figure  5.  The  RMS  differences  for  Cases  1-6  between  the  real  and  the  synthetic  profiles 
computed  using  the  estimated  coefficients  given  in  Fig.  4.  The  dotted  curve  is  the  RMS 
difference  from  Fig.  3b,  shown  for  comparison.  The  shaded  regions  indicate  the  depth  range 
over  which  data  was  deleted. 
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Case  7 


Case  8 


Mode 


Mode 


Case  9 


Case  lO 


Mode 


Mode 


Case  1 1 


Case  12 


Mode 


Mode 


Figure  6.  The  normalized  estimated  coefficient  amplitudes  versus  mode  number  for  Cases 
7  •  12  in  which  data  in  the  lower  portion  of  the  profiles  were  deleted.  Table  3  describes  the 
data  deleted  for  each  case. 
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Case  7 


Case  8 


Case  9 


Case  10 


Case  1 1 


Case  12 
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RMS  Difference  (°C)  RMS  Difference  (°C)  RMS  Difference  (°C) 


Figure  7.  The  RMS  differences  for  Cases  7-12  between  the  real  and  the  synthetic  profiles 
computed  using  the  estimated  coefficients  given  in  Fig.  6.  The  dotted  curve  is  the  RMS 
difference  from  Fig.  3b,  shown  for  comparison.  The  shaded  regions  indicate  the  depth  range 
over  which  data  was  deleted. 
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Case  13  Case  14 


Mode  Mode 


Case  15 


Mode 


Figure  8.  The  normalized  estimated  coefficient  amplitudes  versus  mode  number  for  Cases 
13  -  IS  in  which  data  around  300  m  (a  data  variance  minimum)  were  deleted.  Table  4 
describes  the  data  deleted  for  each  case. 
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Case  13 


Case  15 


Case  14 


RMS  Difference  (°C)  RMS  Difference  (°C)  RMS  Difference  (°C) 


Figure  9.  The  RMS  differences  for  Cases  13  -  IS  between  the  real  and  the  synthetic  profiles 
computed  using  the  estimated  coefficients  given  in  Fig.  8.  The  dotted  curve  is  the  RMS 
difference  from  Fig.  3b,  shown  for  comparison.  The  shaded  regions  indicate  the  depth  range 
over  which  data  was  deleted. 
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Case  16  Case  17 


Mode  Mode 


Case  18 


Mode 


Figure  10.  The  normalized  estimated  coefficient  amplitudes  versus  mode  number  for  Cases 
16  -  18  in  which  data  around  600  m  (a  data  variance  maximum)  were  deleted.  Table  5 
describes  the  data  deleted  for  each  case. 
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Case  16 


Case  17 


Case  18 


RMS  Difference  (°C)  RMS  Difference  (°C)  RMS  Difference  (°C) 


Figure  11.  The  RMS  duierences  for  Cases  16  -  18  between  the  real  and  the  synthetic 
profiles  computed  using  the  estimated  coefficients  given  in  Fig.  10.  The  dotted  curve  is  the 
RMS  difference  fi’om  Fig.  3b,  shown  for  comparison.  The  shaded  regions  indicate  the  depth 
range  over  which  data  was  deleted. 
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Figure  12.  RMS  differences  between  the  data  in  Fig.  la  and  GDEM  climatology  (dotted 
line),  Levitus  climatology  (dashed  Une)  and  synthetic  profiles  computed  from  data  truncated 
below  700  m  to  simulate  T-7  XBTs  (solid  line).  At  all  depths  the  synthetic  profiles  exhibit 
a  smaller  RMS  difference  and  hence  provide  a  better  estimate  of  the  profiles  below  700  m 
than  either  of  the  climatologies. 
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Appendix  A 

FORTRAN  Subprogram 


21 


FORTRAN  Subroutine  ESTCOE 


I.  Introduction 

Subroutine  ESTCOE  estimates  EOF  coeMcients  for  a  set  of  input  data  profiles.  It  also 
calculates  synthetic  profiles  using  the  estimated  coefficients. 

II.  Input 

ESTCOE  requires  the  unit  numbers  of  2  input  files: 
inunitl:  The  mean  profile  &  eigenvectors. 
inunit2:  The  data  profiles  themselves. 

The  mean  profile  &  eigenvectors  must  be  arranged  in  the  first  input  file  as  follows: 

Depth  values  (of  the  mean  profile)  are  stored  in  column  1,  and  the  mean  profile 
values  at  their  respective  depths  must  be  storeu  in  column  2.  An  "end  of  record"  flag 
is  required  at  the  end  of  the  mean  profile.  This  end  of  record  flag  is  a  pair  of  red 
numbers,  9999.00  9999.00. 

After  the  mean  profile  end  record  flag,  the  eigenvectors  are  stored  sequentially. 
For  each  eigenvector,  the  depth  values  should  be  stored  in  column  1,  and  the 
eigenvecto*-  values  at  their  respective  depths  ir.us.  appear  in  column  2.  An  end  of 
record  flag  is  required  at  the  end  of  each  eigenvector  "record." 

NOTE:  The  number  ot  depths  in  the  mean  profile/eigenvector  file  has  been  set  to 
50  via  a  parameter  statement:  parameter(Maxz  =  50).  This  is  also  the  maximum 
number  of  c'V'envectors  allowed  in  the  mean  profile/eigenvector  file. 

The  d:  a  profiles  must  be  arranged  in  the  second  input  file  as  follows: 

Depth  values  are  stored  in  column  1,  and  the  data  values  at  their  corresponding 
depths  are  stored  in  column  2.  Depth  values  should  be  sequential;  that  is,  a 
measurement  taken  at  depth  10  m  should  be  stored  before  a  measurement  taken  at 
20  m.  An  end  of  record  flag  is  required  at  the  end  of  each  profile. 

NOTE:  The  number  of  depths  in  any  single  profile  has  been  set  to  6000  via  a 
parameter  statement:  parameter(Maxpts  =  6000). 

III.  Processing 

First,  the  mean  profile  &  eigenvector  input  file  is  read.  The  user  is  asked  to  decide  how 
many  modes  will  be  used  to  build  the  coefficients  (and  subsequent  synthetic  profiles). 
NOTE:  The  maximum  number  of  modes  available  is  equal  to  the  number  of  eigenvectors 
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in  this  input  file.  Then,  one  at  a  time,  the  profiles  are  read,  coefficients  are  calculated,  and 
synthetic  profiles  are  constructed  using  the  following  approach: 

a.  When  a  profile  is  loaded,  it  is  checked  against  the  eigenvector  depths  for  missing 
data  points  (relative  to  the  eigenvector  depths).  If  a  profile  contains  data  at  depths 
that  are  not  in  the  eigenvector  file,  these  "extra  depffis"  are  ignored  in  subsequent 
coefficient  estimation  subroutines.  If  a  profile  contains  "missing  data"  relative  to  the 
eigenvector  file,  this  is  noted  in  a  "flag"  array. 

b.  For  a  single  profile,  the  system  of  equations  given  in  eq.  (3)  is  solved  using 
standard  Gauss-Jordan  elimination  routines^. 


c.  The  eigenvector  matrix  [e^Cj ...  e^]  is  created  using  only  the  eigenvector  depths 
where  input  data  values  are  present.  To  solve  for  the  coefficient  vector  represented 
by  the  set  of  values  {c„},  the  inverse  of  the  eigenvector  matrix  is  found.  The 
determinant  of  the  resulting  matrix  is  displayed.  If  the  determinant  is  too  small 
(relative  to  a  set  tolerance  of  l.OD-14),  the  user  is  warned  that  there  may  be  large 
errors  in  the  resulting  estimated  coefficients.  The  user  is  encouraged  to  rerun  the 
program  using  a  smaller  number  of  modes  with  which  to  build  the  coefficients. 

d.  The  right-hand  side  vector  in  eq.  (3)  is  created  using  the  input  profile  data,  the 
mean  profile,  and  the  eigenvectors.  The  entire  system  of  equations  is  then  solved. 

e.  The  synthetic  profile  is  built  using  the  estimated  coefficients.  The  coefficients  and 
the  synUietic  profile  are  stored,  and  the  program  loops"  through  again  using  the  next 
proffie’s  data. 

IV.  Output 


Output  from  ESTCOE  consists  of  the  variables  listed  below.  Parameter  values  are  set  for 
the  maximum  number  of  depths  allowed  in  the  mean  profile/eigenvector  input  file  (Maxz 
=  50)  and  the  maximum  number  of  input  profiles  (MaxProfiles  =  100). 


Num__probes: 

e_n_dep: 

Num_eigen: 

e__depth: 


Number  of  profiles  in  the  input  data  file 

Number  of  depths  in  each  eigenvector  (as  well  as  the  mean 

profile) 

Number  of  eigenvectors  (or  modes)  used  in  calculating  the 
estimated  coefficients 

Array  of  eigenvector  depths;  e_depth(MaxDepths) 


^Press,  William  H.  et  al.  (1986).  Numerical  Recipes  -  The  Art  of  Scientific  Computir^. 
New  York  (NY):  Cambridge  University  Press,  pp.  35-37. 
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Coef: 

pp: 


2-D  array  that  stores  the  estimated  coefficients; 
Coef(MaxProfiles,  Maxz) 

2-D  array  that  stores  the  synthetic  profiles;  pp(MaxProfiles, 
Maxz) 
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V. 


Example  of  input  files  and  driver  program 

(1)  Input  file  1:  Mean  temperature  profile  and  eigenvectors 


0.00 

28.38 

(depth  0  m,  mean  profile  value  at  0  m) 

10.00 

28.33 

1750.00 

4.03 

2000.00 

3.66 

(depth  2000  m  ,  mean  profile  value  at  2000  m) 

9999.00 

9999.00 

(end  of  record) 

0.00 

-0.01 

(depth  0  m,  eigenvector  1) 

10.00 

0.04 

1750.00 

0.10 

2000.00 

0.04 

(depth  2000  m,  eigenvector  1) 

9999.00 

9999.00 

(end  of  record) 

0.00 

0.06 

(depth  0  m,  eigenvector  N) 

10.00 

0.04 

1750.00 

-0.19 

2000.00 

0.00 

(depth  2000  m,  eigenvector  N) 

9999.00 

9999.00 

(end  of  record) 

(2)  Input  file  2:  Temperature  profiles 

0.00 

28.60 

(depth  0  m,  temperature  value  at  0  m,  profile  1) 

2.00 

28.60 

4.00 

28.58 

5.00 

28.57 

8.00 

28.40 

10.00 

2824 

••• 

1750.00 

3.96 

9999.00 

9999.00 

(end  of  record) 

b!oo 

28.29 

(depth  0  m,  temperature  value  at  0  m,  profile  N) 

5.00 

28.28 

10.00 

28.27 

1700.00 

3.81 

9999.00 

9999.00 

(end  of  record) 
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(3)  Example  of  a  driver  program  that  calls  the  subroutine  ESTCOE. 


$large 

program  driver 

c  Purpose:  To  test  the  subroutine  ESTCOE 

c  Compiler:  Microsoft  FORTRAN  v  4.01 
c  To  compile:  fl  -c  driver.for 
c  To  link:  link  driver  +  ESTCOE; 
c  To  run:  driver 

c  Parameter  values  set:  Maximum  number  of  depths  in  any  input  file  is  set  at  Maxz 
c  Maximum  number  of  input  profiles  is  set  at  MaxProfiles  =  100 

parameter(Maxz  =  50,  MaxProfiles  =  100) 

double  precision  Coef(MaxProfiles,  Maxz),  pp(MaxProfiles,  Maxz) 

integer  inunl,  inun2,  outunl,  outun2 

integer  Num_eigen,  Num_probes,  e_n_dep,  mode,  probe 

real  e_depth(Maxz) 

c  Set  unit  numbers  for  input/output  files 
inunl  =  10 
inun2  -  11 
outunl  =  12 
outun2  =  13 

c  Give  ii^t  filenames 

open(inunl,  file  =  ’evec.dat’) 
open(inun2,  file  =  ’profles.dat’) 

c  Give  the  output  filenames 
open(outunl,  file  =  ’coeff.dat’) 
open(outun2,  file  -  ’synthetdat’) 

c  Call  the  estimate  coefficient  subroutine 

call  est_coef(munl,  inun2,  Num_probes,  e_n_dep,  Num_eigen, 

•  e_depth,  Coef,  pp) 

c  Write  the  coefficients  to  output  filel 
4  format(lx,  i4,  lx,  flO.4) 

do  10  probe  =  1,  Numj)robes 
do  20  mode  =  1,  Num_eigen 
write(outunl,  4)  mode,  Coef(probe,  mode) 
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20  continue 
10  continue 


c  Write  the  synthetic  profiles  to  output  file2 
5  fomiat(lx,  flO.2,  lx,  flO.4) 
do  30  probe  =  1,  Num_probes 
do  40  k  =  1,  e_n_dep 
write(outun2,  5)  e_depth(k),  pp(probe,  k) 
40  continue 
30  continue 

c  Close  all  open  files 
close(inunl) 
close(inun2) 
close(outunl) 
close(outun2) 

c  End  of  program 
stop 
end 
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$large 

^  mmmmmm****m***m*****m*******m**m****m*m***m***m*m*****m*m*m****m*****m 

subroutine  ESTCOE(inunl,  inun2,  Num_probes,  e_n_dep, 

*  Num_eigen,  e_depth,  Coef,  pp) 

Q  m*mm»**m*******m****’tt**********m*******m*******************mmm***mm**m 


C 

c  Purpose: 
c 
c 
c 

c  Input: 

c 

c 

c  Ou^ut: 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


To  find  EOF  coefficients  for  input  profiles  which 
have  missing  depth  values  and  to  create  synthetic 
profiles. 

inunl:  Unit  number  of  input  profiles 
inun2:  Unit  number  of  Mean  profile  &  Eigenvectors 

Num_probes:  Number  of  profiles  in  the  input  data 
file 

e_n_dep:  Number  of  depths  in  each  eigenvector  (as 

well  as  the  mean  profile) 

Num_eigen:  Number  of  eigenvectors  (or  modes)  used 
in  calculating  the  estimated 
coefficients 

e^depth:  Array  of  eigenvector  depths; 

e_dq)th(Maxz) 

Coef:  2-D  array  which  stores  the  estimated 

coefficients;  Coef(MaxProfiles,  Maxz) 
pp:  2-D  array  which  stores  the  synthetic 

profiles;  pp(MaxProfiles,  Maxz) 


c  Limitations:  Maximum  number  of  input  profiles  is  set  at 
c  MaxProfiles  =  100.  Maximum  number  of  dq)ths  in  any 

c  single  input  profile  is  set  at  Maxpts  =  6(XX). 

c  Maximum  number  of  dq)ths  in  the  mean/eigenvector  file 

c  is  set  at  Maxz  =  SO.  Matrix  Singularity  criteria  set 

c  at  det_tol  =  l.OD-14. 


c 


c  Precisicm:  Double  Precision 


c  Compiler:  Microsoft  FORTRAN  v  4.01 
c  To  compile:  fl  -c  ESTCOE.for 
c  To  link:  link  driver  +  ESTCOE; 


c 


c  Code  uses  routines  fiom  NUMERICAL  RECIPIES 
c  ROUTINES:  L  U  Decomposition  (LUDCMP) 
c  Back  Substitution  (LUBKSB) 

c  NOTE:  Modifications  wav  made  to  the  original  code  for  these  two 
c  subroutines;  real  variables  were  changed  to  double  precision; 
c  unnecessary  ou^ut  variables  were  deleted  fiom  calling  and 
c  subroutine  declmtion  statements. 
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c 

c  Pipgnun  writtra  by:  Eileen  P.  Kennelly 
c  Nq)tune  Sciences,  Inc 

c  SUdeU,  LA  70458 

c  Variables: 

c  A:  "A"  matrix  in  A  *  x  =  B;  Used  in  LUDCMP 

c  B:  "B"  vector  inA*x  =  B;x  =  (inv)A  *  B;  The 

c  coefficients  for  a  particular  profile  are  temporarily 

c  stored  in  this  vector;  Used  in  LUBKSB 

c  Coef:  2-D  array  which  stores  calculated  coefficients; 

c  Coef(MaxProfiles,  Maxz) 

c  determin:  Determinant  of  the  LUDCMP  matrix;  checks  for 
c  singularity 

c  det_tol:  Tolerance  check  of  the  determinant  of  the  LUDCMP 

c  output  matrix;  set  at  parameter  (det_tol  =  l.OD-14) 

c  E:  Eigenvectors  are  stored  in  this  2-D  matrix; 

c  E(Maxz,  Maxz) 

c  eof:  End  of  input  file  flag 

c  e_depth:  Array  of  eigmvector  depths;  e_depth(Maxz) 

c  e_n_dep:  Number  of  depths  in  each  eigenvector 

c  f:  Shallow  water  data  (input  profiles)  is  stored  here 

c  one  profile  at  a  time;  a  1-D  array,  f(Maxpts) 

c  fbar:  Mean  profile  associated  with  the  eigenvectors;  a  1-D 

c  array  fbar(Maxz) 

c  flag:  Logical  array;  ^  to  TRUE  if  there  exists  an 

c  eigenvector  value  at  a  profile  dqith 

c  f^dqnh:  Array  of  profile  dqiths;  f_dq)A(Maxpts) 

c  f_n_dq):  Number  of  dqiths  in  a  particular  profile 

c  fjjrime:  At  dqith  z,  f  =  Profile  value  -  Mean; 

c  f  (z)  *  f(z)  -  fbar(z);  a  ID  array  fj)rime(Maxpts) 

c  indx:  Array  used  by  LUDCMP  and  LUBKSB 

c  inunl:  Unit  number  of  input  Mean  profile  and  eigoivectors 

c  inun2:  Unit  number  of  input  profiles 

c  MaxProfiles:  Maximum  numb^  of  input  profiles; 
c  parameter  (MaxProfiles  =  1(X)) 

c  Maxz:  Maximum  number  of  depths  in  either  input  file; 

c  parameter  (Maxz  =  50) 

c  N:  Number  of  eigenvectors  in  the  eigenvector  file 

c  Num_eigen:  Number  of  eigenvectors  which  are  actually  used  in  the 
c  construction  of  estimated  coefficients 

c  Numjprobes:  Number  of  profiles  in  input  file 
c  p:  1-D  array  which  stores  a  single  synthetic  profile; 

c  p(Maxz) 

c  pp:  2-D  array  which  stores  all  synthetic  profiles; 

c  pp(MaxProfiles,  Max) 

c  probe:  Profile  counter 
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Q  m***m********»*m*mm*******mm*<iHf*******mit‘*m*m***m*mm*m****m*m********** 


double  precision  det_tol 

parameter  (Maxz  =  50,  MaxProfiles  =  100,  Maxpts  =  6000) 
parameter  (det_tol  =  l.OD-14) 

double  precision  A(Maxz,  Maxz),  B(Maxz),  f(Maxpts),  determin 
double  precision  E(Maxz,  Maxz),  fbar(Maxz),  f_prime(Maxpts) 
double  precision  p(Maxz),  Coef(MaxProfiles,  Maxz) 
double  precision  pp(MaxProfiles,  Maxz) 
integer  inunl,  inun2,  indx(Maxz) 

integer  Num_eigen,  probe,  e_n_dep,  f_n_dep,  N,  Num_probes 

logical  eof,  flag(Maxz) 

real  e_dep^(Maxz),  f_depth(Maxpts) 

c  Get  the  Eigenvectors  and  the  Mean  vector  from  first  input  file 
call  EigMea(inunl,  N,  E,  e_n_dq),  e_depth,  fbar) 
call  Ld_Eig(N,  Num_eigen) 

do  10  probe  =  1,  MaxProfiles 

c  Let  user  know  something  is  happening... 

write(*,*)  ’Calculating  coefficients  for  profile  ’,  probe 

c  Load  an  input  profUe,  f 

call  Ld_Pr^inun2,  f,  f_depth,  f_n_dep,  eoO 
if  (eof)  goto  89 

c  Find  missing  depths  and  set  the  flag  array  to  indicate  this 
call  Miss(f_depth,  f_n_dq),  e_depth,  e_n_dep,  flag) 

c  Load  the  matrix  A  (in  A  *  x  «  B)  with  the  Hgravectors 
call  A_Matx(e_n_dq),  E,  Num^eigoi,  flag.  A) 

c  Call  the  Numerical  Recipe  LUDCMP  to  find  the  inverse  of  the  A 

c  matrix  in  A  *  x  =  B 

call  LUDCMP(Num_eigai,  A,  indx) 

c  Look  at  the  determinant  of  the  invose  A  matrix.  Check  it  for 

c  singularity 

call  Det(Num_eigen,  A,  determin) 
if  (dabsfdeterinin)  .It.  d^_tol)  thoi 
call  MessageO 
endif 

c  Calculate  fjirime;  f  (z)  =  f(z)  -  fbar(z) 

call  Ld_lpr(f,  f_d^th,  f_n_dqp,  e_dq)th,  e_n_dq). 
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*  fbar,  f_priine) 

c  Load  vector  B  (in  A  *  x  =  B)  with  f  *  E 

call  Ld_B(E,  e_depth,  e_n_dep,  f_depth,  f_n_dep, 

•  fj)rime,  Num_eigen,  B) 

c  Solve  for  x:  A  *  x  =  B 
c  (resulting  coefficients  stored  in  variable  B) 
call  LUBKSB(A,  Num_eigen,  indx,  B) 

c  Create  the  synthetic  profile, 

caU  Syn(B,  E,  fbar,  e_n_dep,  Nuni_eigen,  p) 

c  Store  coefficients  in  a  2-D  array 

caU  Mk_Coe(probe,  Num_eigen,  B,  Coef) 

c  Store  synthetic  profile  in  a  2-D  array 
caU  Mk_Syn(probe,  e_n_dep,  p,  pp) 

10  continue 

89  Nuni_probes  =  probe  -  1 
return 
end 

subroutine  EigMeaGun,  N,  E,  e_n_dep,  e_dq>th,  fbar) 

c  Purpose:  To  load  the  Mean  vector,  and  the  Eigenvectors 

parameter  (Maxz  »  SO) 

double  precision  E(Maxz,  Maxz),  fbar(Maxz) 

int^er  lun,  N,  tjijiep 

real  e_dq)th(Maxz) 

c  Load  the  Mean  vector,  fbar 
c  End  of  mean  profile  flag  set  at  9999.00 

do  10  k  =  1,  Maxz 
read(lun,  *)  e_dq)th(k),  fbar(k) 
if  (e_dq>th(k)  .gt.  9000.  .and.  fbar(k)  .gt.  9000.)  goto  20 
10  continue 
20  e_n_dq)  =  k  -  1 

c  Load  the  eigenvectors,  E(eigenmode,  depth) 
c  End  of  each  eigmvector  flag  set  at  9999.00 
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do  30  j  =  1,  Maxz 
do  40  k  =  1,  Maxz 

read(lun,  *,  end  =  50)  e_depth(k),  EO,  k) 
if  (e_dq)thOc)  .gt.  9000.  .and.  k)  .gt.  9000.)  goto  31 
40  continue 

31  continue 

30  continue 

c  Note  that  there  are  (j  -  1)  eigenvectors  available 
50  N  =  j  -  1 

return 

end 

subroutine  Ld_Eig(N,  num_eigen) 

c  Purpose:  The  user  decides  how  many  modes  will  be  used  to  build 
c  the  coefficients  (and  eventually  the  synthetic  profiles). 

parameter  (Maxz  =  50) 
integer  N,  num_-ig^.n 

10  format(lx,  A,  13,  lx,  A,  A) 

20  write(*,  10)  ’There  are’,  N,  ’eigenvalues  in  the  ’, 

*  ’input  eigenvector  file.  You  can  use  ’ 
write(*,  10)  ’a  maximum  of,  N,  ’eigenvalues  in  the  ’, 

*  ’construction  of  the  synthetic  profile.  ’ 
write(*,*)  ’How  many  do  you  want  to  use?  ’ 
read(*,*)  num_eigen 

return 

end 

subroutine  Ijd_Pro(lun,  f,  f_dq)th,  f_n_dep,  eof) 

c  Purpose:  To  load  a  single  profile 

parameter  (Maxz  =  50,  Maxpts  =  6(XX)) 

double  precision  f(Maxpts),  vals(2,  Maxpts) 

integer  lun,  f_n_dq) 

logical  eof 

real  f_depth(Maxpts) 

c  Load  a  proffle,  f 
c  End  of  profile  flag  set  at  9999.00 
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do  10  k  =  1,  Maxpts 

read  (lun,*,  end  =  40)  vals(l,  k),  vals(2,  k) 
if  (vals(l,  k)  .gt.  9000.  .and. 

*  vals(2,  k)  .gt.  9000.)  goto  20 
f_depth^)  =  vals(l,  k) 
f(k)  =  vads(2,  k) 

10  continue 

20  f_n_dep  =  k  -  1 
goto  30 

40  eof  =  .true. 

30  continue 

return 

end 

subroutine  Miss(f_depth,  f_n_dep,  e_depth,  e_n_dep,  flag) 

c  Purpose;  To  find  the  missing  points  in  a  profile 

parameter  (Maxz  =  SO,  Maxpts  =  6000) 
integer e_n_dq),  f.n_dep 
logical  fiagCMaxz) 

real  e_depth(Maxz),  f_depth(Maxpts) 

do  10  k  =  1,  e_n_dep 
flag(k)  =  .false, 
do  20  nn  =  1,  f_n_dep 
if  (f_dq)th(nn)  .eq.  e_dq)th(k))  then 
flag(k)  =  .true, 
goto  10 
endif 

20  continue 
10  continue 

return 

end 

^  m*mm*ti‘****m**m*mm******m***********************************m********** 

subroutine  A_Matx(e_n_dep,  E,  Num_eigen,  flag,  A) 
c  Purpose:  To  load  the  "A"  Matrix  in  A  •  x  =  B 
parameter  (Maxz  =  50) 

double  precision  E(Maxz,  Maxz),  A(Maxz,  Maxz) 
int^er  e_n_dep,  Num_eigai 
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logical  flag(Maxz) 
real  sum 


c  Create  the  A  matrix  using  only  the  eigenvectors  needed, 
do  10  i  =  1,  Num_eigen 
do  20  j  =  1,  Num_eigen 
sum  =  0.0 

do  30  k  =  1,  e_n_dep 
if  (flag(k))  then 

sum  =  sum  +  E(i,  k)  *  E(j,  k) 
endif 

30  continue 

A(i,  j)  =  sum 
20  continue 
10  continue 

return 

end 

subroutine  Det(N,  A,  prod) 

c  Purpose:  Check  the  determinant  of  the  (inv)A  matrix  for 
c  Dossible  singularity 

parameter  (Maxz  =  50) 

double  precision  A(Maxz,  Maxz),  prod 

integer  N 

prod  =  1.0 
do  10  i  =  1,  N 
prod  =  prod  ♦  A(i,  i) 

10  continue 

4  format(lX,  A,  IX,  E8.2) 

write(*,4)  ’Determinant  of  the  Eigenvector  LUDCMP  matrix  =’, 
*  prod 


return 

Old 

^  ti******)**********************^^**************************************** 

subroutine  Ld_fipr(f,  f_depth,  f_n_dep,  e_depth,  e_n_dep, 

*  fbar,  fjprime) 

c  Purpose:  At  each  depth  z,  calculate  f  (z)  =  f(z)  -  fbar(z) 
parameter  (Maxz  =  50,  Maxpts  =  6000) 


double  precision  fbar(Maxz),  f_prime(Maxpts),  f(Maxpts) 

integer  f_n_dep,  e_n_dep 

real  e_depth(Max2),  f_depth(Maxpts) 

do  10  k  =  1,  e_n_dep 
do  20  nn  =  1,  f_n_dep 
if  (f_depth(nn^.eq.  e_depth(k))  then 
f_prime(nn)  =  f(nn)  -  fbar(k) 
endif 

20  continue 
10  continue 

return 

end 

Q  ^c^cicitintmmmm^m^enirti***********************************************^!******** 

subroutine  Ld_B(E,  e_depth,  e_n_dep,  f  depth,  f_n_dep, 

*  f_prime,  Num_eigen,  B) 

c  Purpose;  Load  the  vector  "B"  in  A  *  x  =  B 

parameter  (Maxz  =  50,  Maxpts  =  6000) 

double  precision  B(Max2),  f_prime(Maxpts),  E(Maxz,  Maxz),  sum 

integer  Num_eigen,  f_n_dep,  e_n_dqp 

real  e_dq)th(Maxz),  TdepthCMaxpts) 

do  10  i  =  1,  Num_eigen 
sum  =  0.0 

do  20  k  =  1,  e_n_dq) 
do  30  nn  =  1,  f_n_dep 
if  (f_depth(nn)  .eq.  e_depth(k))  then 
sum  =  sum  +  f_prime(nn)  ♦  E(i,  k) 
endif 

30  continue 
20  continue 

B(i)  =  sum 
10  continue 

return 

end 

subroutine  Syn(C,  E,  fbar,  e_n_dep,  Num_eigen,  p) 

c  Purpose:  At  each  depth  z,  calculate  the  synthetic  profile,  p 
c  p(z)  =  Mean(z)  +  sum  (Coefficient  *  Evec) 

c  Note:  This  is  for  a  single  proftle 
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parameter  (Maxz  =  50) 

double  precision  C(Maxz),  fbar(Maxz),  E(Maxz,  Maxz),  p(Maxz),  sum 
integer  Num_eigen,  e_n_dep 

c  Synthetic(z)  =  Mean(z)  +  sum  (Coeff  *  Eigenvector) 
do  20  k  =  1,  e_n_dep 
sum  =  0.0 

do  10  j  =  1,  Num_eigen 
sum  =  sum  +  C(j)  *  EO,  k) 

10  continue 

p(k)  =  fbar(k)  +  sum 
20  continue 

return 

end 

subroutine  Mk_Coe(probe,  Num_eigen,  B,  Coef) 

c  Purpose:  To  store  the  calculated  coefficients  in  a  2-D  matrix 

parameter  (Maxz  =  50,  MaxProfiles  =  100) 
double  precision  B(Maxz),  Coef(MaxProfiles,  Maxz) 
integer  probe,  Num^eigen 

do  10  k  =  1,  num_eigen 
Coef(probe,  k)  =  B(k) 

10  continue 

return 

end 

Q  m**mm*n‘*****m***********m********************************************* 
subroutine  Mk_Syn(probe,  e_n_dep,  p,  pp) 

c  Purpose:  To  store  the  calculated  synthetic  profiles  in  a  2-D  matrix 

parameter  (Maxz  =  50,  MaxProfiles  =  1(X)) 
double  precision  p(Maxz),  pp(MaxProfiles,  Maxz) 
integer  probe,  e_n_dep 

do  10  k  =  1,  e_n_dep 
pp(probe,  k)  =  p(k) 

10  continue 

return 

end 

Q  mm***0m*****m******mmm***#******m*mtif**m*»****m**mm***m*#***********m* 
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subroutine  LUDCMP(N,  A,  indx) 


c  Purpose;  To  solve  a  Matrix  A  into  L  U  Decompostion  form 
c  ...  from  NUMERICAL  RECIPES 

parameter  (Maxz  =  SO) 

parameter  (nmax  =  100,  TINY  =  l.OE-20) 

double  precision  A(Maxz,  Maxz),  VV(nmax),  aamax,  sum,  dum 

integer  N,  indx(Maxz),  D,  imax 

do  26  i  =  1,  N 
indx(i)  =  0 
26  continue 


D=1 

do  12i=l,N 
aamax =0. 
do  11  j  =  l,N 

if  (dabs(A(i,  j))  .gt.  aamax)  aamax  =  dabs(A(i,  j)) 

1 1  continue 

if  (aamax  .eq.  0.)  pause  ’Singular  matrix.’ 

W(i)  =  1.  /  aamax 

12  continue 

do  19  j  =  1,  N 
if  (j  .gt.  1)  then 
do  14i  =  l,j  -  1 
sum  =  A(i,  j) 
if  (i  .gt.  1)  then 
do  13  k  =  1,  i  -  1 
sum  =  sum  -  A(i,  k)  *  A(k,  j) 

13  continue 
A(i,j)  =  sum 

endif 

14  continue 
endif 

aamax  =  0. 
do  16  i  =  j,  N 
sum  =  A(i,  j) 
if  (j  .gt.  1)  then 
do  15  k  =  1,  j  -  1 
sum  =  sum  -  A(i,  k)  *  A(k,  j) 

15  continue 
A(i,  j)  =  sum 

endif 

dum  =  W(i)  *  dabs(sum) 
if  (dum  .ge.  aamax)  then 
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imax  =  i 
aamax  =  dum 
endif 

16  continue 

if  (j  .ne.  imax)  then 
do  17  k  =  1,  N 
dum  =  A(imax,  k) 

A(imax,  k)  =  A(j,  k) 

A(j,  k)  =  dum 

17  continue 
D=-D 

W(imax)  =  W(j) 
endif 

indx(j)  =  imax 
if  0  .ne.  N)  then 

if(A0,j)  .eq.  0.)  AO,  j)  =  TINY 
dum  =  1.  /  AO,  j) 
do  18  i  =  j  +  1,  N 
A(i,  j)  =  A(i,  j)  *  dum 

18  continue 
endif 

19  continue 

if  (A(N,  N)  .eq.  0.)  A(N,  N)  =  TINY 

return 

end 

subroutine  LUBKSB(A,  N,  indx,  B) 

c  Purpose:  To  solve  the  system  A  *  x  =  B  for  x  by  back  substitution 
c  ...  from  NUMERICAL  RECIPES 

parameter  (Maxz  =  SO) 

double  precision  A(Maxz,  Maxz),  B(Maxz),  sum 
integer  indx(Maxz),  11,  N 

ii  =  0 

do  12  i  =  1,  N 
11  =  indx(i) 
sum  =  B^) 

B(U)  =  B(i) 
if  (ii  .ne.  0)  then 
do  1 1  j  =  ii,  i  -  1 
sum  =  sum  -  A(i,  j)  *  BO) 

1 1  continue 

elseif  (sum  .ne.  0.)  then 
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ii  =  i 
endif 

3(i)  =  sum 

12  continue 

do  14  i  =  N,  1,  -1 
sum  =  B(i) 
if  (i  .It.  N)  then 
do  13  j  =  i  +  1,  N 
sum  =  sum  -  A(i,  j)  *  B(j) 

13  continue 
endif 

B(i)  =  sum  /  A(i,  i) 

14  continue 

return 

end 

subroutine  MessageQ 

c  Purpose:  To  flash  a  warning  message  if  the  determinant  of  the 
c  LUDCMP  matrix  is  too  small 

write(*,*)  ’  ’ 

write(*,*)  ’WARNING.  In  trying  to  calculate  coefficients,  it 

*  ’was  found  that  the’ 

write(*,*)  ’eigenvector  matrix  was  singular;  the  coefficients 

*  ’may  be  erroneous.  ’ 

write(*,*)  ’It  is  suggested  that  you  re-run  this  program  and  ’, 

*  ’choose  a  smaller  number’ 

write(*,*)  ’of  eigenvalues  to  construct  the  coefficients.’ 
write(*,*)  ’Processing  continues...’ 
write(*,*)  ’  ’ 


return 

end 
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